{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 260,
   "metadata": {
    "tags": [
     "remove-cell"
    ]
   },
   "outputs": [],
   "source": [
    "# Reference: https://jupyterbook.org/interactive/hiding.html\n",
    "# Use {hide, remove}-{input, output, cell} tags to hiding content\n",
    "\n",
    "import sys\n",
    "import os\n",
    "if not any(path.endswith('textbook') for path in sys.path):\n",
    "    sys.path.append(os.path.abspath('../../..'))\n",
    "from textbook_utils import *"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Minimizing Huber Loss\n",
    "\n",
    "*Huber loss* combines absolute loss and squared loss to get a function that is differentiable (like squared loss) and less sensitive to outliers (like absolute loss):\n",
    "\n",
    "$$\n",
    "L(\\theta, \\textbf{y}) = \\frac{1}{n} \\sum_{i=1}^n \\begin{cases}\n",
    "    \\frac{1}{2}(y_i - \\theta)^2 &  | y_i - \\theta | \\le \\gamma \\\\\n",
    "     \\gamma (|y_i - \\theta| - \\frac{1}{2} \\gamma ) & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "Since Huber loss is differentiable, we can use gradient descent. We first find the gradient of the average Huber loss:\n",
    "\n",
    "$$\n",
    "\\nabla_{\\theta} L(\\theta, \\textbf{y}) = \\frac{1}{n} \\sum_{i=1}^n \\begin{cases}\n",
    "    -(y_i - \\theta) &  | y_i - \\theta | \\le \\gamma \\\\\n",
    "    - \\gamma \\cdot \\text{sign} (y_i - \\theta) & \\text{otherwise}\n",
    "\\end{cases}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We create the functions `huber_loss` and `grad_huber_loss` to compute the average loss and its gradient. We write these functions to have signatures that enable us to specify the parameter as well as the observed data that we average over and the transition point of the loss function:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 258,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-01T02:28:40.909807Z",
     "start_time": "2018-07-01T02:28:40.905060Z"
    }
   },
   "outputs": [],
   "source": [
    "def huber_loss(theta, dataset, gamma=1):\n",
    "    d = np.abs(theta - dataset)\n",
    "    return np.mean(\n",
    "        np.where(d <= gamma,\n",
    "                 (theta - dataset)**2 / 2.0,\n",
    "                 gamma * (d - gamma / 2.0))\n",
    "    )\n",
    "\n",
    "def grad_huber_loss(theta, dataset, gamma=1):\n",
    "    d = np.abs(theta - dataset)\n",
    "    return np.mean(\n",
    "        np.where(d <= gamma,\n",
    "                 -(dataset - theta),\n",
    "                 -gamma * np.sign(dataset - theta))\n",
    "    )"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we write a simple implementation of gradient descent. The signature of our function includes the loss function, its gradient, and the data to average over. We also supply the learning rate. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 259,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-01T02:28:40.869020Z",
     "start_time": "2018-07-01T02:28:40.864546Z"
    }
   },
   "outputs": [],
   "source": [
    "def minimize(loss_fn, grad_loss_fn, dataset, alpha=0.2, progress=False):\n",
    "    '''\n",
    "    Uses gradient descent to minimize loss_fn. Returns the minimizing value of\n",
    "    theta_hat once theta_hat changes less than 0.001 between iterations.\n",
    "    '''\n",
    "    theta = 0\n",
    "    while True:\n",
    "        if progress:\n",
    "            print(f'theta: {theta:.2f} | loss: {loss_fn(theta, dataset):.3f}')\n",
    "        gradient = grad_loss_fn(theta, dataset)\n",
    "        new_theta = theta - alpha * gradient\n",
    "        \n",
    "        if abs(new_theta - theta) < 0.001:\n",
    "            return new_theta\n",
    "        \n",
    "        theta = new_theta"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall that the bus delays dataset consists of over 1,000 measurements of how many minutes late the northbound C-line buses are in arriving at the stop at 3rd Avenue and Pike Street in Seattle:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 261,
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "delays = pd.read_csv('data/seattle_bus_times_NC.csv')"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In {numref}`Chapter %s <ch:modeling>`, we fit a constant model to these data for absolute loss and squared loss. We found that absolute loss yielded the median and square the mean of the data:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 262,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Mean:   1.920\n",
      "Median: 0.742\n"
     ]
    }
   ],
   "source": [
    "print(f\"Mean:   {np.mean(delays['minutes_late']):.3f}\")\n",
    "print(f\"Median: {np.median(delays['minutes_late']):.3f}\")    "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we use the gradient descent algorithm to find the minimizing constant model for Huber loss: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 264,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2018-07-01T02:28:41.054144Z",
     "start_time": "2018-07-01T02:28:40.912188Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Minimizing theta: 0.701\n",
      "\n",
      "CPU times: user 93 ms, sys: 4.24 ms, total: 97.3 ms\n",
      "Wall time: 140 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "theta_hat = minimize(huber_loss, grad_huber_loss, delays['minutes_late'])\n",
    "print(f'Minimizing theta: {theta_hat:.3f}')\n",
    "print()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The optimizing constant for Huber loss is close to the value that minimizes absolute loss. This comes from the shape of the Huber loss function. It is linear in the tails and so is not affected by outliers like with absolute loss and unlike with squared loss."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ":::{warning}\n",
    "\n",
    "We wrote our `minimize` function to demonstrate the idea behind the algorithm. In practice, you will want to use well-tested, numerically sound implementations of an optimization algorithm. For example, the `scipy` package has a `minimize` method that we can use to find the minimizer of average loss, and we don't even need to compute the gradient. This algorithm is likely to be much faster than any one that we might write. In fact, we used it in {numref}`Chapter %s <ch:donkey>` when we created our own asymmetric modification of quadratic loss for the special case where we wanted the loss to be greater for errors on one side of the minimum than the other.  \n",
    "\n",
    ":::"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "More generally, we typically stop the algorithm when $ \\theta^{(t)} $ doesn't\n",
    "change much between iterations. In our function, we stop when $ \\theta^{(t+1)} - \\theta^{(t)} $ is less than 0.001. It is also common to stop the search after\n",
    "a large number of steps, such as 1,000. If the algorithm has not arrived at the\n",
    "minimizing value after 1,000 iterations, then the algorithm might be diverging\n",
    "because the learning rate is too large or the minimum might exist in the limit\n",
    "at $ \\pm \\infty $. \n",
    "\n",
    "Gradient descent gives us a general way to minimize average loss when we cannot\n",
    "easily solve for the minimizing value analytically or when the minimization is\n",
    "computationally expensive. The algorithm relies on two important properties of\n",
    "the average loss function: that it is both convex and differentiable in $ \\boldsymbol{\\theta} $. We discuss how the algorithm relies on these properties\n",
    "next."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.4"
  },
  "toc": {
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
